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Abstract: We present a Monte-Carlo (MC) calculation of the propagation of cosmic ray protons in the 
Galaxy for energies above 1 PeV. We discuss the relative strengths of competing effects such as paral- 
lel/perpendicular diffusion and drifts in toy models of the Galaxy. We compare our estimates with the 
results of the MC calculation for the toy models and then we apply the MC calculation to a few more 
realistic models of the Galactic magnetic field. We study the containment times in different models of the 
magnetic field in order to understand which one may be consistent with the low energy data. 
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Introduction 

The transport of cosmic rays in the Galaxy is usu- 
ally studied solving the diffusion-convection equa- 
tion from a distribution of sources in a medium 
with given diffusion properties. This approach, al- 
though very successful (see Ref . [ 1 ] for a recent re- 
view), can not be applied at arbitrarily high energy 
because at some point the diffusion approximation 
breaks down. In the Galaxy this happens, for pro- 
tons, around 10^^ eV. The study of CR transport 
in this energy region is performed using numerical 
simulations of the propagation of single particles 
in the Galactic magnetic field (GMF). Besides the 
obvious advantage of being able to study the transi- 
tion region, another advantage of the latter method 
is the ability to use more realistic models of the 
GMF, including the arms, all the gradients and so 
on. The big drawback of this approach is in the 
computing time that makes it usable only at high 
energy (10^^ eV and above). 

From observations at low energy, ^ GeV, the resi- 
dence time of particles in the Galaxy is estimated at 
about lO*^ yr with an energy dependence of E^'^'^. 
When extrapolated to higher energies, this trend 
produces nonsensical results, predicting, for ex- 
ample, huge anisotropics already around 10^^ eV 
where they are not observed. Assuming an en- 
ergy dependence of E^^^"^, as one would expect 
from Kolmogorov turbulence, the extrapolations 



are much better, but still a few times larger than 
the experimental results [2] . 

In principle one would want to use the diffusion- 
convection equation at low energy, the numerical 
simulations at high energy and possibly match the 
two results in between. Up to now, the simulations 
at high energy [3] were not successful in this re- 
spect. The containment time at 10^^ eV is calcu- 
lated as 10^ yr with an energy dependence of 1/E. 
Extrapolated to lower energy, this result produces 
escape times exceeding the age of the universe. 

In order to investigate this discrepancy, we devel- 
oped a numerical simulation of the propagation of 
particles in arbitrary magnetic fields and we ap- 
plied it to several toy models of the Galaxy to sep- 
arately study the various effects contributing to the 
transport, such as parallel and perpendicular diffu- 
sion, drifts and so on. For details on the simulation 
code see Ref. [4] . 

Diffusion Coefficients 

We calculate the diffusion coefficients in a mag- 
netic field composed of a constant regular field in 
the z direction and a turbulent field with a Kol- 
mogorov spectrum. We inject particles in this field 
and we follow their trajectories recording their po- 
sitions as a function of time. Fitting the distribution 
of the particle positions as a function of time we 
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calculate the diffusion coefficients along the three 
axes. For the parallel diffusion coefficient we find 
the usual energy dependence of E'^/'^ at low energy 
and E"^ at high energy, whereas for the perpendicu- 
lar one we find a steeper dependence at low energy: 
^0.5-0.6 -pjjjg jg already an interesting result be- 
cause it means that in particular geometries where 
the transport is dominated by perpendicular diffu- 
sion one can obtain, from a Kolmogorov spectrum 
of turbulence, a diffusion coefficient, and therefore 
a residence time, steeper than the usual one. 

Toy Models 

The first toy model we consider has a purely az- 
imuthal magnetic field that is constant, 1 fxG, ev- 
erywhere. Superimposed to this regular field there 
is a turbulent component whose magnitude is pro- 
portional to the one of the regular field with a pro- 
portionality constant of 0.5, 1 or 2. We inject parti- 
cles at the solar system position (8.5 kpc from the 
center) and we record the time required for escap- 
ing a cylinder with half-height of 0.5 kpc and ra- 
dius of 10 kpc. It is important to note that in this 
scenario the transport is dominated by the perpen- 
dicular diffusion and possibly by drifts, and that 
the parallel diffusion can be completely neglected 
since it scatters the particles back and forth along 
the field lines, but it does not help them escap- 
ing the disk since the field lines are closed. The 
times of escape are plotted with colored points in 
Fig. 1 . The red upward, blue downward triangles 
and green squares correspond to three different lev- 
els of turbulence: SB/B = 0.5, 1, 2 respectively'. 
The thick solid straight line represents the drift 
time-scale calculated for the average drift veloc- 
ity, whereas the thin solid lines are the time-scales 
for perpendicular diffusion; they are proportional 
to h^/D±. The perpendicular diffusion coefficient 
is the one calculated as described in the previous 
section. As it is clear from Fig. 1 the transport de- 
scription as perpendicular diffusion is pretty good 
for the two cases SB/B — 1 and 2, but it is not 
very accurate for SB /B = 0.5. The discrepancy 
is at low energy, around lO'^^ eV, and in this re- 
gion the drift time-scale is much bigger than the 
diffusion one and then this effect is not caused sim- 
ply by drifts: they are relevant only at higher ener- 
gies. The effect is however produced by the cur- 
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Figure 1: Times of escape in the azimuthal field. 
The colored points are the results of the simula- 
tions for different scenarios. The thin lines are ex- 
pected escape times calculated with the diffusion 
coefficient, whereas the thick straight line is the 
timescale for drifts. 

vature of the field lines since injecting the particles 
at 85 kpc, where the curvature is ten times smaller, 
produces a time of escape that is in much better 
accord with the diffusion time-scale (see light blue 
upward triangles in Fig. 1). 

In order to better understand this effect we calcu- 
late the diffusion coefficients in the curved mag- 
netic field to check if the curvature of the field 
lines has an effect also on the diffusion coefficients. 
Our results are plotted in Fig. 2. We find that with 
smaller levels of turbulence the two perpendicular 
diffusion coefficients, in the radial and in the z di- 
rection, are modified. The first one is increased and 
the second one is decreased. This is consistent with 
the results of Fig. 1, since a smaller diffusion coef- 
ficient in the z direction produces a larger time of 
escape. Reducing the curvature of the field lines, 
i.e. calculating the diffusion coefficients at larger 
distances from the center, the effect is reduced (see 
dotted lines in Fig. 2). The effect promptly disap- 
pears increasing the level of turbulence. Fig. 2 is 
then telling us that the curvature of the field lines 
not only produces a drift, or the rigid displacement 
of the distribution of the particle positions, but it 
also affects the way the distribution evolves with 
time along the three directions. 

Though keeping in mind that this is only a toy 
model of the magnetic field of the Galaxy, it is in- 

1 . Please note that here and in the following when de- 
noting SB we actually mean: {SB^). 
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Figure 2: Diffusion coefficients in a curved reg- 
ular field for the level of turbulence indicated in 
the plots. In each plot the black lines are the par- 
allel (upper) and perpendicular (lower) diffusion 
coefficients obtained in the constant field. The 
green, red and blue lines are diffusion coefficients 
respectively in the azimuthal, radial and z direc- 
tion. Solid lines are for injection at 8.5 kpc, dotted 
lines are for injection at 85 kpc. 

teresting to note that at energy 10^^ eV the escape 
time is T15 w 0.5 — 5 million years (the halo height 
here is only 0.5 kpc). These numbers are of the 
same order of magnitude of the confinement times 
estimated from the abundance of light element in 
the GeV region, which means that in order to fit 
these observations one should postulate that the es- 
cape time below 10^^ cV should be practically en- 
ergy independent. We could not envision any real- 
istic mechanism able to justify such an expectation. 
It follows that within the limitations of the present 
toy model it is very hard to obtain a realistic, even 
qualitative, description of what is observed in the 
Galaxy at much lower energies. This conclusion is 




Figure 3: Galactic Magnetic Field and drifts in the 
BSS models. The black dot is the galactic center, 
the red one is the solar system position. The black 
lines show the position of the center of the arms. 
The blue and red arrows represent the magnitude 
and direction of the GMF in the given position, 
while the black arrows show the drift speed. 

confirmed also by the curves on the grammage that 
show that at 10^^ eV cosmic rays traverse already 
a column density of 1 — 2 g cm~^. 

We modified the above toy model introducing sev- 
eral complications mimicking the structure of the 
Galaxy, for example a gradient in radial direction 
or in the z direction, but the general result did not 
change much. The slopes of the times of escape 
remained quite steep, although the absolute values 
became smaller due to the increase of the drifts. 

"Realistic" Galactic Magnetic Fields 

The knowledge of the magnetic field in our Galaxy 
is very poor, both for what concerns its regular 
component and more so for what concerns the tur- 
bulent one. Having this in mind we considered a 
very simple model with the intent not of being very 
realistic, but just as an example of how the simple 
picture presented in the previous section becomes 
suddenly much more complicated. Indeed, follow- 
ing Ref. [5, 6], we considered the BSS model: it 
includes the spiral arms, has a radial gradient and 
a gradient in the z direction, but it does not have 
an halo field, a bar in the center or other complica- 
tions. 

We plot the drift velocity calculated in several 
points along the line connecting the galactic cen- 
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Figure 4: Times of escape in the BSS model for 
two different ways of normalizing the turbulent 
component of the field. See text for details. 



ter (black dot) with the solar system position (red 
dot) in Fig. 3. The black lines show the position of 
the center of the arms, the blue and red arrows rep- 
resent the magnitude and direction of the GMF in 
the given position, while the black arrows show the 
drift speed. In the toy model presented in the pre- 
vious section the drift velocity was always in the z 
direction and had a position dependent magnitude. 
In this case the situation is much more complicated 
with the drift pushing the particles sometimes to- 
ward the center of the arms and sometimes toward 
the inter-arms space. Moreover in this case both 
the magnitude and direction of the drift depend on 
the position and on the pitch angle of the particle. 
Indeed, in the plot of Fig. 3 the drift velocity is cal- 
culated for a particle with an injection cosine with 
respect to the local magnetic field of \/2/3 and for 
a position just above the galactic plane. Calculat- 
ing it for another pitch angle or another position 
would produce completely different results. For 
example, considering the same pitch angle, but a 
position just below the plane, we would obtain the 
opposite sign for the radial component of the drift 
velocity and thus a very different picture. 



In Fig. 4 we plot the residence times of protons in- 
jected at the solar system position and collected at 
a cylinder with an half-height of 4 kpc and a radius 
of 20 kpc. In each panel we plot several series of 
points corresponding to different levels of turbu- 
lence as indicated. The two panels correspond to 
two different normalizations of the turbulent com- 
ponent of the field. In the upper panel the (SB^) 
of the turbulent field is in every position a fraction 
of the regular field. This produces an almost neg- 
ligible field in the space between the arms, where 
the regular field switches direction. In the lower 
plot the magnitude of the turbulent field follows 
the radial and z dependence of the regular field, 
except for the arms that are absent in the turbulent 
component. This means that in this case the ratio 
5B /B is not constant, but variable, being the value 
indicated in the plot in the center of the arms and 
going to infinity in the space between the arms. 

The times of escape reflect the increased complex- 
ity of the model: they are much less regular than 
the ones obtained in the toy model. The normal- 
ization is however of the same order of magnitude 
with times of escape at 10^^ cV of about 10^ yr and 
grammages of a few gcm^^. The slopes for the 
times of escape are ~ 0.7 and there is no sign of a 
possible transition to a flatter slope that would help 
reconcile these values with the ones measured at 
low energies. 
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